砂丘とかでみられる波波のやつ、風紋の数理モデルがあるらしいので簡単に実装してみた。
風紋モデルの運動には3つのステップがあるようで、saltation, suspension, creep という感じらしい。
saltationは日本語で跳躍とか跳動とか訳されるようで、砂粒が風で飛ばされてしばらくして着陸するみたいな運動。今回のモデルではランダムな点で決まった回数だけ跳躍させた。
suspension（浮遊）は黄砂みたいなフワフワと長く飛ばされる現象で、扱いが面倒なのでモデルには入れない。
creepはそのまんまクリープで、風やら重力やらで表面を転がるのを表している。このモデルでは高さの拡散効果みたいなもんだと思う。
二次元の高度の場を $h(x,y)$とする。
風は$X$正の方向に吹くとして、砂の飛距離を、$L=L_0+bh(x,y)$ で表す。$L_0$は風の効果、$bh$は高さの効果（高いところの砂粒のほうが飛ぶでしょということ）。
あとは簡単で、飛ばされた砂の量(高さ)を $q$とすると、saltation が起こったセル(砂が運ばれる前と後のセル)では
$$ h'(x,y)=h(x,y)-q, \quad h'(x+L(h(x,y)),y)=h(x+L(h(x,y)),y)+q $$と高度が変化する。$h'$は、次の時間ステップに行く前の中間ステップ。
拡散方程式を有限差分法で解くみたいなもん。省略。
上のふたつを交互に発展させて計算するぞ！初期値は適当にランダムノイズを乗せとく。
境界条件は周期境界、計算領域は100$\times$100で、時間発展は1000ステップにしておいた。
using Plots
function main()
dh=0.3
nx=100
ny=100
q =0.1
b =1.0
L0=2.0
D =0.1
h=dh*rand(ny,nx).+2.70 # initial condition
hnew=copy(h)
N=1000 # timestep
freq=10
Hs = Array{Array{Float64,2}}(undef, N÷freq+1)
num=1
Hs[num]=copy(h)
for itr=1:N
# saltation
for k=1:nx*ny
ir=rand(1:nx)
jr=rand(1:ny)
if h[jr,ir]>0
L=Int(floor(L0+b*h[jr,ir]+0.5))
ib=ir+L
h[jr,ir]-=q
h[jr,ifelse(ib>nx, ib%nx, ib)]+=q
end
end
# creep
for i=1:nx
for j=1:ny
hN=h[ifelse(j==ny,1,j+1),i]
hS=h[ifelse(j==1,ny,j-1),i]
hE=h[j,ifelse(i==nx,1,i+1)]
hW=h[j,ifelse(i==1,nx,i-1)]
hNE=h[ifelse(j==ny,1,j+1),ifelse(i==nx,1,i+1)]
hSE=h[ifelse(j==1,ny,j-1),ifelse(i==nx,1,i+1)]
hSW=h[ifelse(j==1,ny,j-1),ifelse(i==nx,1,i+1)]
hNW=h[ifelse(j==ny,1,j+1),ifelse(i==1,nx,i-1)]
hnew[j,i]=(1-D)*h[j,i]+D*((hN+hS+hE+hW)/6+(hNE+hSE+hSW+hNW)/12)
end
end
# for i=1:nx
# for j=1:ny
# h[j,i]=hnew[j,i]
# end
# end
@.(h=copy(hnew))
if itr%freq==0
num+=1
Hs[num] = copy(h)
print(itr, " ")
end
end
return Hs
end
@time H=main();
10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 960 970 980 990 1000 0.406244 seconds (33.68 k allocations: 9.676 MiB, 5.38% compilation time)
N=length(H)
anim = @animate for t=1:N
h=heatmap(H[t],clim=(0,10),c=:sandyterrain,title="height")
plot(h, size=(500,450))
end
filename="fuumon.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\****\****\****\****\fuumon.gif
うおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおおお
なんか波ができとりますわ。きれい。
N=length(H)
anim = @animate for t=1:N
h=surface(H[t], zlims=(0,30),legend=:none,c=:summer)
plot(h, size=(800,800))
end
filename="fuumon3D.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\****\****\****\****\fuumon3D.gif
